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Abstract 



We compare different nonlinear approximations to gravitational clustering in the weakly 
nonlinear regime, using as a comparative statistic the evolution of non-Gaussianity which 

can be characterised by a set of numbers Sp describing connected moments of the density 
field at the lowest order in (5^): (5"')c — Generalizing earlier work by Bernardeau 

(1992) we develop an ansatz to evaluate all Sp in a given approximation by means of a 
generating function which can be shown to satisfy the equations of motion of a homogeneous 
spherical density enhancement in that approximation. On the basis of the values of we show 
that approximations formulated in Lagrangian space (such as the Zeldovich approximation 
and its extensions) are considerably more accurate than those formulated in Eulerian space 
such as the Frozen Flow and Linear Potential approximations. In particular we find that 
the nth order Lagrangian perturbation approximation correctly reproduces the first n + 1 
parameters We also evaluate the density probability distribution function for the different 
approximations in the quasi-linear regime and compare our results with an exact analytic 
treatment in the case of the Zeldovich approximation. 

Subject Headings: cosmology: theory - large scale structure of Universe 
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1. Introduction 



Gravitational instability in the Universe can be characterised by two distinct epochs: during 
the first, density fluctuations 5 = {p — Po)/Po evolve self-similarly according to the tenets of 
linear theory ((5''^^(x, t) oc D+(t)A(x)) with the result that a density distribution that was 
initially Gaussian remains Gaussian at subsequent epochs as long as \6\ << 1. The linear 
epoch clearly cannot continue indefinitely since a stage will arise when S, although small, 
becomes comparable to unity, so that weakly nonlinear effects can no longer be ignored. 
This is the quasi-hnear regime; as long as |5| < 1 the evolution of the density field can be 
adequately described by means of a perturbative expansion: S{x,t) = E~=i5^"Hx,t) where 
5*^"^(x, t) = D"(t)A*^"-'(x) if the Universe is spatially fiat and matter dominated (Peebles 
1980, Fry 1984). This epoch witnesses the growth of skewness and kurtosis and other higher 
moments of the one-point probability distribution function (PDF) of density perturbations 
characterising a non-Gaussian development of an initially Gaussian field S^^\ During still 
later stages, 5 becomes larger than unity, with the result that perturbative approximations 
break down and a fully nonlinear treatment of the problem is required. This later epoch is 
characterised by the formation of pancakes (5 — > oo along two dimensional sheets) and the 
gradual development of cellular structure (Shandarin & Zeldovich 1989). 

Although no exact treatment is available which encompasses the entire nonlinear epoch, 
some approximations have been suggested which attempt to mimick certain features of non- 
linear gravitational instability. Our aim in this paper is to investigate the relative accuracy 
of five such approximations in the weakly nonlinear regime {\S\ < 1) where they may be an- 
alytically compared to the exact solution in the form of a perturbative expansion in powers 
of 6. 

Our treatment overlaps with (and considerably extends) the recent work of Munshi & 
Starobinsky (1993) (hereafter MS) and Bernardeau ct al. (1993), in which 3 approximations: 
the Zeldovich approximation (Zeldovich 1970), hereafter ZA, the frozen fiow approximation 
(Matarrese et al. 1992), hereafter FF, and the linear potential approximation (Brainerd et 
al. 1993; Bagla & Padmanabhan 1994), hereafter LP, were compared to the exact solution in 
second and third order in perturbation theory. In this paper, we provide an ansatz whereby 
nonlinear approximations can be compared with the exact perturbative solution to any order 
in perturbation theory. We apply this ansatz to approximations formulated in Eulerian space 
such as FF and LP, as well as to approximations based on Lagrangian perturbation theory (of 
which ZA happens to be the lowest order solution) . It appears that Lagrangian perturbative 
methods are, as a rule, much more accurate than either FF or LP to any given order in 
perturbation theory. 

2. Nonlinear approximations 

Gravitational instability in a spatially fiat matter dominated FRW Universe before caustic 
formation can, in the Newtonian approximation, be described by the following system of 
equations (see, e.g., Zeldovich & Novikov 1983, Peebles 1980): 

= AnGa'^Pod; (la) 



3 



S + -W^. ((l + 5)u) = 0; (16) 
a 

(au)- + (uVx)u = -Vx$ (Ic) 

where u = ax, x being the comoving coordinate (a(t) oc t^^^ and po = l/QnGt"^). If we 
assume that u is irrotational then we can define a velocity potential such that u(x, t) — 
Vx'y(x,t) = V[/(x, t) where wc have introduced the notation V = ^Vx, H = a/a = 2/3t 
(the potential U is related to the potential V used in MS by the relation U = —HV). 
Moreover, At/ = ^ div^xXi = 9 which is the dimensionless velocity divergence used in MS. 
In this notation, we can rewrite the equations as follows: 

d 

a— 5(x, a) + (1 + (5(x, a))^(x, a) + V(5(x, a) VC/(x, a) = 0; (2a) 
oa 

d 1 

a— ^(x, a) + -e(x, a) + VC/(x, a) Ve(x, a)) + UikUik + A$(x, a) = 0; (26) 

A$=^5(x,a) (2c) 

where Uik = {-^Y&^ikU and summation over repeated indexes is assumed. The initial 
condition for the system (2) at i — > is given by the linear approximation corresponding to 
a growing adiabatic mode: 

$ = 0o(x); C/ = -^</.o(x); 5=^A0o(x); 0^-8. (3) 

Several of the nonlinear approximations which we will consider such as the Zeldovich 
approximation, the Frozen Flow approximation and the Linear Potential approximation arise 
as a result of the extrapolation of one of the linearised relations in (3) into the nonlinear 
regime. Thus in FF, the value of the velocity potential is kept fixed to its linearised value 
U = — |0o, so that u = — |V0o- The movement of a particle in FF is such that the particle 
upgrades its velocity with each time step to that prescribed by the local value of the linear 
velocity field. This results in a laminar fiow for the velocity field since different particle 
trajectories can never intersect in principle. 

The linear potential approximation on the other hand, is based upon the assumption 
that the gravitational potential does not evolve with time so that $ = 0o(x). As a result 
particles effectively move along the lines of force of the primordial potential 0o- Similarly the 
Zeldovich approximation is based on extending [/ = — 1$ into the nonlinear regime. In all 
three cases the constraint equations U = — 1$ (ZA), U — — |0o (FF), ^ = 0o (LP) replace 
the Poisson equation (2c) which is not satisfied in these approximations beyond the linear 
regime. 

In nonlinear approximations formulated in Lagrangian space, the main object of study is 
the particle trajectory. In these approximations the initial comoving (Lagrangian) coordinate 
q and the Eulerian coordinate of a particle x(t) are related by a displacement field VE^: 

x = q+*(t,q). (4) 
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If we introduce the matrix 

Qq.k dqk 

then Mjfe satisfies the equations 

eikiMkruM-j = 0, (66) 

where J = |det[Mjfe]| is the Jacobian of the transformation (4), eiki is the unit totally antisym- 
metric tensor, Eq. (6b) being the condition of potential motion in Eulcrian space (see, e.g., 
Zeldovich & Novikov 1983). Density and velocity fields are determined from the relations 

^ = 6 + l^J-\ u = a|^] =a*. (7) 

A solution of Eqs. (6a,b) in the quasi-linear regime may be obtained by expanding in 
a power series = "i/^^^ + 'i/^'^^ + .. where higher orders in ^^"^ (n > 1) are related to lower 
orders via an iterative procedure (Moutarde et al. 1991, Buchert 1992, Lachieze-Rey 1993). 
It appears that for the matter dominated Universe all factorize: 

= Dl{t)4''\ci), D^it) = ^ cx a{t). (8) 

Subsequent Lagrangian approximations then arise if one truncates this scries after a finite 
number of terms. As a result of this truncation Eq. (6a) becomes approximate and the 
Poisson equation (la) ceases to be valid, but Eqs. (7) as well as the continuity equation (lb) 
continue to be exactly satisfied in these approximations. 

The first-order Lagrangian approximation ^ = is simply the Zeldovich approxima- 
tion: 

x. = <i. + D^it)4'\<l), #) = -^. (9) 

The account of the next, second order terms results in what one may call the post- Zeldovich 
approximation (hereafter PZA). It is given by * = -|- where 

<? = -nK<')'-^.'M'). (10") 
= </',!■?, (106) 

coma means partial derivative with respect to q and summation over repeated indexes is 
assumed (Bouchet et al. 1992, Gramann 1993 and others). In the third order [post-post- 
Zeldovich approximation, hereafter PPZA), * = Z^n=i ^^"^ with 
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V'g'-€' = ^(^S<>-^S<.') (116) 

(cf. Juszkiewicz et al. 1993, Bernardeau 1993) and so on. A distinguishing feature of PPZA 
is that motion becomes non-potential in Lagrangian space (but not in Eulerian space, of 
course) beginning from this approximation. 

Our aim in the present paper will be to determine how well the five nonlinear approxima- 
tions discussed above perform when compared within the framework of perturbation theory 
which is well defined in the quasi-linear regime when \S\ < 1. To investigate statistical be- 
haviour in all orders of perturbation theory we follow Bernardeau (1992) who developed an 
elegant ansatz (which we summarise below) by means of which vertice weights characterising 
irreducible moments of 6 may be evaluated at every order by means of a generating function 
Gs. 

3. Generating function for the irreducible moments of S. 

Let us introduce a vertice generating function for any random field F(x, a) as 

GF(r) = ^ ii-^r- (12a) 

n=l ^■ 

where 

, _ /(FW(x, a)(5(^)(xi, a)...^(^)(x„, a)),d'^d^^^..A^ 

^ (/((5«(x,a)5(i)(x',a))d3x#x')" ' ^ ^ 

F^^^ is the n-th order of expansion of F in a power series with respect to 5, 5*-^-' is the 
linear approximation for b given in Eq. (3), and only connected diagrams are taken into 
account, which explains the subscript c. Throughout the present analysis we assume that the 
initial density field 5^^) and the associated linear gravitational potential 0o(x) are Gaussian 
stochastic quantities though the results obtained below may be generalized to a non-Gaussian 
case as well. Higher vertices of b and the dimensionless velocity divergence 9 = /\U are 
denoted by z/„ and /i„ respectively so that 

G8 = f: -,r\ Ge = f: ^r^ (13) 

(i/„ = (5*^"))c, Hn = (^*^"^)c with ui — 1^ — —1). We define r with the opposite sign as 
compared to Bernardeau (1992) to simplify the appearance of some expressions. 

It is well known that for small values of a"^ = {6^^^'^) the connected moments of the density 
field have the simple form (Fry 1984, Bernardeau 1992): 

(Oc - SniSy-' (14a) 

where Sn are related to the vertice weights z/„ introduced earlier by 

5-4 = 4i/3 + 12i/| 
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Gp^pH^^ir) = -Gafag- (15o?) 



5*5 = 5i/4 + 60i/3i/2 + 60i/| 

Sq = 6i/5 + 120i/4i^2 + 90i/^ + 720i/3i/^ + 360i/^ (146) 

The equations relating T„ and have the same form as (14b) with T„ replacing Sn and 
being replaced by (— 1)"'^„. The additional (—1)"" factor arises because yU„ are defined with 
the use of 5^^^ according to Eq. (12b). Formulas for the moments of 9 would be completely 
identical to those for 5 if we defined /x^, using 9^^^ = —6^^\ This just corresponds to the 
multiplication of /x„ by (— 1)". 

Now we demonstrate how the vertice weights Vn (and hence Sn) can be determined in any 
arbitrary order for the nonlinear approximations described above. We do this by following 
Bernardeau (1992), who studied leading (tree level) diagrams in the limit when the density 
variance cr^ was very small and showed that the generating functions (12a) defined for any 
arbitrary stochastic fields F and H possess the following properties: 

d 

G^apir) = t—Gf, (15a) 

da (jj- 

Gfh{t) = Gf{t)Gh{t), (156) 
GvF-VH = 0, (15c) 

1 

-( 

3 

Using (15a - d) we can rewrite equations (2a,b) in terms of vertice generating functions 
of the corresponding fields (as a result, they become equations for statistically averaged 
quantities) : 

T-^Gs + {l + Gs)Ge^O, (16a) 

ar 

7 -I 1 

T— + -Gg + -Gl + Ga^ = 0. (166) 

The Poisson equation (2c) is now replaced by: 

(1) Ga$ = for the exact solution, 

(2) Ga$ = —2GAU — ~2^o for ZA, 

(3) Gg = -T for FF, 

(4) Ga$ = Ga0o = 1^ for LP- 

It is not straightforward to write the corresponding conditions replacing the Poisson equa- 
tion for PZA, PPZA and higher Lagrangian approximations because they are formulated 
in Lagrangian space, and Eqs. (16a,b) originated from Eulerian equations. Remarkably, it 
appears unnecessary because as we shall show in the next section it is possible to circumvent 
this problem entirely. 

4. Pcirticle trajectories in the spherical model and values of the moments Umpt-n- 

For the exact perturbative solution, the system of equations (16a,b) may be transformed 
into a second order differential equation for the variable y — 1 + Gs{r). For ZA, FF and 
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LP, Eq. (16b) may be solved directly after substituting the corresponding expressions for 
Ga<s>- However, a much simpler and more physical method exists that is also applicable in 
the cases of PZA, PPZA and higher Lagrangian approximations. 

Let us take a closer look at Eqs. (16a,b). It is easy to see that if we make the substitution 
Gs S, Go ^ 9 and t ^ a, then the statistical Eqs. (16a,b) become dynamical equations 
describing the isotropic and homogeneous expansion of the Universe (a) under the influence 
of gravity in the case of the exact solution, and (b) due to some other forces mimicking 
gravity for the different approximations (because the Poisson equation is not satisfied in this 
case). The initial conditions Gs — t, Gg — —r as r ^ can be satisfied simultaneously if 
a free scaling dimensional coefficient of proportionality in dynamical solutions is chosen so 
that 5 = a for a — > 0. Therefore, the dynamics in question is actually that of the spherical 
top- hat model. This agrees well with an intuitive picture that for a ^ 1, large and rare 
fluctuations {a <^ 5 < 1) are approximately spherical. The same evidently refers to the 
cases of PZA, PPZA and so on, inspite of more comphcated forms for the third equation 
replacing the Poisson equation, because it follows from Eqs. (15c, d) that to obtain statistical 
equations for generating functions one has to neglect inhomogeneities in dynamical equations 
and substitute all second-rank tensors by those proportional to a unit tensor. 

Now, the easiest way to solve the top-hat model is to work with Lagrangian (comoving) 
coordinates q. Then the problem reduces to one of finding appropriate particle trajectories 
x(q, a) associated with equations (16a,b) in each of the approximations and in the exact 
solution respectively. 

The exact spherical top-hat solution: 

Let r = |x|, ro = |q| = r(0). The gravitational potential inside a homogeneous sphere is 
$ = ^9^. It clearly follows from mass conservation that 1 -|- 5 = (ro/r)^. So the equation of 
motion of a particle (a spherical shell) is 



(the independent variable is changed from t to a). The solution of Eq. (18) in a parametric 
form is 




(17) 



This is none other than the usual Newtonian equation = —-^r^ for the reduced physical 
scale R — ar/rQ. Its first integral satisfying the initial condition 5 = a for a ^ is 




(18) 
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(19a) 
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(196) 



2 (1 - cosOf 



2 




(19c) 
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(cf. Zeldovich & Novikov 1983, Peebles 1980). The corresponding expression for Ge follows 
from Eq. (16a). 

The vertice weights Vn are found by expanding Gs{t) near r = 0. As a result, we obtain: 
= T + 0.810r2 + O.eOlr^ + 0.426r^ + 0.293r^ + (20a) 

Ge = -(r + 0.619t^ + 0.376t^ + 0.226t^ + 0.135t^ + ...). (206) 

Comparing (20a) with (13) we find: 1/2 = 34/21 ~ 1.619, i/3 = 682/189 ~ 3.608, i/4 ^ 10.23. 
Higher moments of the density distribution such as the skewness 5*3, the kurtosis 84^ etc. can 
now be determined by substituting the values of into Eq. (14b). 

The Zeldovich Approximation: 

In this case, particle trajectories are already given by Eq. (9). The initial gravitational 
potential 0o(^o) satisfying the initial condition 5 = a for a — > is 0o = a^r^/Qt^ (it doesn't 
depend on t since a{t) oc t^l^). As a result 

r(a)=ro(l-^), (21a) 

^z.^C^.= (^)''-l=(l-0"^-l, (216) 

G, = -t(i-0"\ (21c) 
Expanding G5 and Ge near r = 0, we get 

Gs = T + 0.667t^ + 0.370t^ + 0.185t^ + 0.086t^ + (22a) 

Ge = -(r + 0.333r^ + O.Ulr^ + 0.037r^ + 0.012r^ + ...), (226) 

from which we can readily obtain the values of z/„, /i„ and Sn, Tn- For n = 3, 4, they coincide 
with those found previously (Grinstein & Wise 1987, MS 1993, Bernardeau et al. 1993). 

PZA, PPZA and higher order Lagrangian approximations: 

Here the Lagrangian method shows its superiority over the Eulerian one because particle 
trajectories are explicitly given by Eqs. (4, 8-10) in the case of PZA, Eqs. (4, 8-11) for PPZA 
and so on. Just as in the previous paragraph, the initial gravitational potential for the top- 
hat spherical model is 0o(ci) = a^rg/9t^ where Tq = |q|. Then ip^^^ = — 2a^/3t^ = —a/D^{t), 
L)2 ^(2) = -aVo/21 (the index r means a radial component). Therefore, a particle trajectory 
in PZA is 

, , / a a^\ , , 

r{a) = ^0 (^1 - 3 - , (23a) 

and the generating functions for PZA are 
SpzA = ^5^(l~^~^) -1 = T + 0.810t2 + 0.561r^ + 0.358r^ + 0.215t^ + (236) 



Ge = = -(r + 0.619t2 + 0.254t3 + 0.114r^ + O.OSOr^ + ...)■ (23c) 

In this case, the first two coefficients of the expansion (20a,b) are reproduced exactly as a 
result of which PZA gives the correct value for the skewness in the lowest order in a. 

Furthermore, D\'il)f> = -23a^ro/1701 (see Eq. (11a)), so that a PPZA particle trajectory 
is given by 

/ a 0,2 23a'^\ , 

ria) = ro 1 . 24a 

^ ' \ 3 21 1701 j ^ ' 

Thus, the generating functions for PPZA have the forms: 

/ 2 23r^\"^ 

G^= 1 — - — — - 1 = r + O.SlOr^ + O.eOlr^ + 0.412t^ + 0.268r^ + (246) 

= \^ ^^3^3 = -(r + 0.619t^ + 0.376t^ + 0.168t^ + 0.082t^ + ...). (24c) 

~ 3 ~ 21 ~ 1701 

Now the first three terms in the expansion (20a,b) arc correctly reproduced. Therefore, the 
PPZA values for both the skewness 6*3 as well as the kurtosis 5*4 are exact. 

It is clear now that in the A^-th order Lagrangian approximation defined by the condition 
* = S^=i a particle trajectory in the spherical top-hat model has the form r(a) = 
To (1 — Pn{ci)) where is a A^-th order polynomial with Pn{0) — 0, so the generating 
functions have the following structure: 

Gs^il-PN{r))-'-l, (25a) 

Ge ^ -3^^ {1 - Pj.{T)r\ (256) 
ar 

The expansion of Gs, Gq around the point r = correctly reproduces the first N terms of 
the series for the exact solution (20a,b) and consequently the first A" + 1 values of Sn- 

The Frozen Flow Approximation 

This is the only case when it is easier to solve the problem in Eulerian space because the 
velocity potential is equal to its initial, Gaussian distributed value. Therefore, 9 is Gaussian, 
too, and Gq = — r. Consequently 

Gs = expij) - 1. (26) 

However, the calculation using a particle trajectory method is not long either. The 
velocity potential in the top-hat model in FFA is equal io U — —2(j)o{r)/3 — —2a^r'^/21t^. 
So, the equation of motion of a particle has the form 

Ur = a—^ -2a'^r/9t. (27a) 
dt 
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Its solution in terms of a is (Matarrese et al. 1992) 

r(a) = ro exp(-^) (276) 

that just gives the expression (26) for the density generating function. 

The Linear Potential Approximation 

In LP, it is the gravitational potential $ that stays equal to its initial linear value 0o(r) 
which in our case is once more equal to a^r^/9t^ . The equation of motion of a particle is 
therefore 

d ( r,dr\ 2a^ 



(cf. Eq. (17)) or, it terms of a. 



da^ 2a da 2a 
Solving Eq. (28b) we obtain (cf. Brainerd et al. 1993) 



a^-l ) = -Zr_r (28a) 
dt\ dt 9^2 ^ ^ 



d'^r 3 dr r „ 

33 + 7r:3r + 7r: = 0. (286) 



sin V2a. (29a) 




- 1 = T + 0.567t^ + 0.242t2 + 0.087t^ + 0.028t^ + (296) 

\s'iny zr / 

G0 = ^ (v^coiv^ -i)^ -(r + 0.133t2 + 0.025t^ + O.OOSlr^ + O.OOlOr^ + ...) (29c) 
in LP. 

In figure 1 we plot the density contrast in the spherical top-hat model obtained using 
PPZA, PZA, ZA, LP and FF against the exact solution. We find that for a positive density 
perturbation all models underestimate the density contrast, but PPZA, PZA and ZA are 
more accurate than LP and FF. This is quite surprising since it is well known that the 
Zeldovich approximations is the least accurate during spherical collapse when all eigenvalues 
of the deformation tensor d'^U/dqidgj are equal (Shandarin & Zeldovich 1989). We also list 
the turnaround epoch and the recollapse epoch (in units of r) in each of the approximations 
and in the exact solution in Table 1. We find that both turnaround and recollapse occur 
later in the approximations than in the exact solution. 

Expanding and Gq in each of the approximations near r = we get the reduced 
moments z/„ yU„ and, from (14b), the parameters Sn and T„. Our results are summarised in 
Table 2 for the first six moments of the distribution: 5*1... 5*6. Table 3 contains the values of 
Ti...Tq. The values which we obtain for 5*3 and 5*4 for ZA, FF and LP are identical to those 
obtained earlier by MS and Bernardeau et al. (1993). Interestingly we find that of all the 
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approximations considered by us, PPZA appears to be the most accurate in reproducing the 
results of the exact perturbative treatment. Next in accuracy comes PZA which reproduces 
the skewness only, then ZA followed by LP and FF. The results of our analysis lead us to 
conclude that nonlinear approximations formulated in Lagrangian space (ZA, PZA etc.) are 
significantly more accurate in the weakly nonlinear regime than those which are formulated 
in Eulerian space (FF and LP). It would be of great interest to extend this analysis to the 
strongly nonlinear regime where the different approximations should be compared with re- 
sults of the adhesion model and N-body simulations. Some work in this direction is presently 
in progress (Sathyaprakash et al. 1994). 

5. Density distribution functions 

Having obtained the generating functions Gg, it is fairly straightforward to calculate cor- 
responding density PDFs ri{5) for each of the approximations in the limit o" <^ 1 and for 
a sufficiently small \5\ (for the exact perturbative solution this was done by Bernardeau 
(1992)). A more detailed discussion of the region of applicability of the approach used in 
the paper will be given below. 

Let V{y) be the Laplace transform of ri{5) (with a displayed explicitly): 

V{y) = j\{5)exp (-^^^) dS . (30) 

Then the generating function of the moments Sn (see Eq. (14a)): 

m = T.Sp-^y''^ ^1 = ^2 = 1, (31) 
p=i p- 

is connected to V by the simple formula 'P{y) = exp{—Lp{y)/a'^) (see, e.g.. White 1979). On 
the other hand, as can be checked by a direct comparison of the coefficients of the series 
(13) and (31) using (14b) that (j){y) can be expressed through the function C(r) = 1 + Gs{t) 
using the following relations: 

^{y)=yCir{y)) + ^T\ (32a) 

r{y)^-yC{r{y)) (326) 

where dot means derivative with respect to r (see, e.g., Bernardeau & Schaeffer 1992), so 
that 

f0i-c%-c^(r-fy\ (320 

Finally, the density PDF follows from the inverse Laplace transformation 




(72 



dy (33) 



where c = const. 
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For example, for a purely Gaussian stochastic field, ({t) = 1 + t, y = —r and ipiy) = 
y — 7/^/2 in accordance with the fact that all Sp with p > 2 in Eq. (31) are equal to zero in 
this case. Then the inverse Laplace transform of this ip{y) is just the Gaussian distribution 
r]{5) = (27r(72)-V2exp (-572(72). 

In the limit a <^ 1, the main contribution to the integral in (33) is made by stationary 
points of the exponent in it, i.e. by roots of the equation 

^^l + G,(r) = l + 5 (34) 
dy 

(complex roots of this equation should be considered, too). Thus, here we return to the 
prescription 5 ^ Gs used in the derivation of equations for the generating functions in Sec. 
3. In this (steepest descent) approximation, Eq. (33) takes the form 




where now t{5) should be determined from Eq. (34) for any given generating function Gs{t). 

Let us now discuss the accuracy of the formulae (33, 35) for the density PDF. Since 
the generating functions Gs{r) and ^p{y) have been calculated in the lowest order in cr 
only, there exist small corrections to them beginning from terms proportional to so that 
ip{y) = (po{y) + cr'^^iiy) + ••• (here the zero-order term (pQ{y) stands for ip{y) used above) 
and similarly for Gs- The account of the correction results in the multiplication of 
the integrand of Eq. (33) by the cr-independent term exp{—(fi{y)) having an unknown 
dependence on y. Similarly, the multiplicative factor exp{—Lpi{y{T))) will appear in the 
right-hand side of Eq. (35). In the limit \6\ <C 1 when 5 pa r ~ —y, this factor has the 
form (1 + 0(5^)), the 5^ (or r^) correction arising due to a cr^ correction to the dispersion 
of density perturbations. Of course, we may renormalize by defining it to be equal to the 
exact dispersion of density perturbations and not {S^^^'^) as is used in the paper, but all other 
corrections beginning from a term oc will remain. 

Therefore, we have to conclude that if each term in the power series (13) for the generating 
function Gs is computed in the lowest order in cr as is done in Bernardeau 1992 and in 
the present paper, then the resulting density PDF has an exponential accuracy: the large 
exponential factor in it appears to be correct for not too large d, but the coefficient of the 
exponential should be taken in the limit 5—^0. One is allowed to keep terms linear in 6 in the 
latter but we shall not do so except for retaining the coefficient I/VQjw^ in front of the large 
exponential (we make an exception for the Zeldovich approximation where it is instructive to 
keep the entire coefficient in Eq. (35) in order to compare it with an exact analytic solution). 
This is enough, however, to obtain very significant deviations from Gaussian behaviour in 
the quasi-linear regime. 

On the other hand, to go beyond exponential accuracy and evaluate the coefficient of 
the exponential exactly, one has to calculate (pi{y). Generally, this quantity is spectrum- 
dependent although the ZA presents an important exception to this rule. But even then, the 
resulting expression may not be used for arbitrarily large 5 in the limit a <^ 1. In particular. 
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it is not possible to get the threshold value Sc appearing in a Press-Schechter-like formula for 
the number of collapsed objects from such a treatment. This arises because other density 
configurations different from the spherical top-hat one may give a dominant contribution to 
the generating function Gs in this regime. It results for instance, in that the value 5c for ZA 
should be taken as a/5 in the limit cr <^ 1 as shown in the Appendix, and not as 3 as seems 
from Eq. (38) below. Moreover, the case of the Zeldovich approximation shows that the 
value of 6 for which such a new dominant configuration appears cannot be obtained from 
the spherical approximation for the vertice generating function. 
Now we consider the concrete approximations. 

The Zeldovich approximation: 

This case is very important because we may compare our result with an exact analytic 
solution for the density PDF valid for all values of a and 5 which was obtained by Kofman 
(1991). This permits us to check the range of validity of the method used. This solution and 
its limit for o" <^ 1 are given in the Appendix. 

Using Eqs. (21b, 32a,b, 33) and changing the integration variable in (33) from y to r, 
we obtain 

({r)=(l-^)\ y = -T(l-^)\ ip{y) = -T + lT\ (36a) 



3/ ' " V 3/ ' 6 



where the integration contour C corresponding to that in Eq. (33) is a continuous curve in 
the complex r plane beginning in the sector |r| oo, O.Ttt < argir) < O.Qtt and ending in 
the sector |t| — > oo, — 0.97r < arg{T) < — O.Ttt. The only irregular non-analytic point of 
the integrand is at |r| = cxo. So, for cr <^ 1, only stationary points of the exponent should 
be considered. There are four of them: 



A0n 



2n{n - 1) 



""- = 3(^1- (jqr^J' ^=1,2,3, /3n = ^3 n = Q.Q. (37) 

The last point does not make a contribution in the lowest order in a because of the coefficient 
in front of the exponential in Eq. (36b). The point Ti gives the largest exponent for not too 
large S. The direction of the steepest descent for it is perpendicular to the real axis. Thus, 
we find using (35) that 



,W = ^(l-5) (l--) e-i^,r(,)^3ll-j—j-\. (38) 



Comparing this expression with the rigorous result for ZA in the hmit a <^ 1 (Eq. (A. 2)), 
we see that the large exponential is reproduced exactly up to a rather large value of S but 
the coefficient of the exponential in Eq. (38) is correct in the limit of small \S\ only. In 
the latter limit, Eq. (38) just reduces to Eq. (A. 3) if terms up to 0{6) inclusive are kept 
in the coefficient. 0(6'^) and higher terms in series expansions of the coefficients of the 
exponential in (38) and (A. 2) are different. Moreover, the appearance of the square root 
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singularity in Eq. (38) at r = 0.6, 6 ~ 0.95 is an artifact of the approach used, this point 
is completely regular in the rigorous density PDF for ZA. All this is in complete agreement 
with the general discussion of limits of validity of the spherical approximation for the vertice 
generation function in the beginning of this section. (Note also that the exponent in Eq. 
(38) coincides with that found by Padmanabhan & Subramanian (1993) for the PDF of the 
final smoothed density field in ZA using a completely different approach.) 

PZA and PPZA: 

With the accuracy chosen above, we immediately get from Eqs. (23b, 24b, 34, 35) that 



where 




(39a) 



for PZA and 

^ir)=(l-l-"^-^^] -1 (396) 



-2 OQ-^3\-3 



3 21 1701 



for PPZA. 



The Frozen flow approximation: 

Now it follows from Eqs. (26, 34, 35) that r = ln{l + 5) and 

where we have introduced the coeficient (1 + 6)^^ (that is possible within our accuracy) to 
make the PDF being normalized to 1 exactly. PDF (40) is just the log-normal distribution 
which was independently proposed as a good statistical approximation for the density PDF 
in the quasi-linear regime by Hamilton (1988) and Coles & Jones (1991) (in contrast to the 
dynamzcfl/ approximations which we are considering). The result (40) shows that there exists 
a close internal relationship between FFA and the log-normal approximation. In particular, 
we may expect that they have approximately the same accuracy in the quasi-linear regime. 

The Linear Potential approximation: 

Eqs. (29b, 34, 35) are now relevant, giving the result 

6. Conclusions 

The nonlinear evolution of an initially Gaussian distribution leads to the following relation- 
ship between the connected density moments (5^"'^)c at the lowest order in (5^): (5")c — 
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Sn{S'^)^~^ where the parameters Sp characterise the development of non-Gaussianity (5*3 de- 
scribes the skewness of the density distribution and 5*4 its kurtosis). Generahsing earher 
work by Bernardeau (1992) we have shown that for approximation methods attempting to 
mimick the effects of nonhnear gravity the values of Sp can be derived by means of a generat- 
ing function Gs- For a given nonlinear approximation Gs = 5sph is simply the overdensity in 
the spherical top hat model in that approximation. Thus knowing how particle trajectories 
evolve in the spherical top hat model in a given nonlinear approximation we can determine 
Gs and consequently Sp and the probability distribution function for the density. 

Following this ansatz we determine the functional form of the generating function and 
the values of the first six parameters 5*1, ...Sq for five distinct nonlinear approximations, of 
which three are formulated in Lagrangian space (including the Zeldovich approximations and 
its extensions), and the remaining two in Eulcrian space (frozen flow and linear potential). 
Comparing our results with those of an exact perturbative treatment (Bernardeau 1992) we 
find that an approximation formulated to n*'^ order in Lagrangian space correctly reproduces 
the first n -\- 1 parameters Sn (see Table 2). Our comparison leads us to conclude that 
nonlinear approximations which are formulated in Lagrangian space are considerably more 
accurate than those formulated in Eulerian space when tested in the weakly nonlinear regime. 
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APPENDIX 



The rigorous formula for the density PDF in the Zeldovich approximation for a single-stream 
motion in Eulerian space with Gaussian initial conditions is (Kofman 1991, see also Kofman 
et al. 1994) 

"^^^ = A.aHl + 5)^ ' + e-i^+e ^ - e ^ ds, (Ala) 



47ra\l + 5)3 Jso{5) 
3 

:i+^ 



So{S) = .^fM/.v Pn{s,S) =s-V5^ + COS ^{n- 1)71+ ^arccos^(^^^ - ij^j . 

{A.lb) 

For (J <S 1 and not too large 6, the main contribution to the integral comes from the lower 
limit of integration s — Sq. If s = so{l + x), x <^ 1, then Pi — Sq ■ 3-\/5/2 + 0{x), P2,3 — 
+so ■ \/hx (1 =F + 0{x)). Then the terms with /?2,3 are the leading ones (they almost 

cancel each other), and the formula (A. la) takes the form 

If \5\ <C 1, too (but \5\ may be much more than cr), then (A. 2) simplifies to 

The expression (A. 2) clearly looses sense for 5 > (7/2)3 _ 1 ~ 42, but it is no more 
the dominant term in (A. la) for this value and even for somewhat smaller values of S. 
A numerical calculation shows that for 6 > 30.6, the main contribution to the PDF is 
produced by the maximum of the exponential exp (— ((s — 3)^ + /5|) /2a'^) which is located 
at s ~ 4/3 for 5 — > 00, other terms in (A. la) being exponentially smaller. As a result, 
r){S) oc (1 + 5)~^exp{-5^/2a'^) with 5c = a/S for 5 ^ 00 and a <^ 1. Thus, the Press- 
Schechter-like density perturbation value Sc (a "threshold" for formation of compact objects) 
is equal to a/5 in the Zeldovich approximation, and not to 3 that would follow from a 
naive application of the formula (A. 2) in the regime 5 3> 1. Of course, even \/5 is a great 
overestimation of a real value of Sc in the exact solution (if it has sense at all) which is 
expected to be equal to, or a little bit less than 1.686. 

Note also that the PDF (A. la) is not normalized to unit total probability due to the 
appearance of multistreaming motion in the Zeldovich approximation (as well as in the 
exact solution) even for a small a. Actually, J!^iri{6)d6 gives the mean number of streams 
Ng (see also Kofman et al. 1994). This is, however, an exponentially small effect, it is 
straightforward to show by direct integration of Eq. (A. la) that 

for (7 <^ 1. Here, the threshold value Sc — appears once more. This formula is not bad 
even for cr = 1 where it gives Ng ~ 1.0087 instead of the exact value A^s(l) ~ 1.0137. 
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Exact 


PPZA 


PZA 


ZA 


LP 


FF 




1.062 


1.118 


1.194 


1.5 


2.058 


3 


Tcoll 


1.686 


2.050 


2.266 


3. 


4.935 


oo 



Table 1: Spherical collapse 





Exact 


PPZA 


PZA 


ZA 


LP 


FF 




4.857 


4.857 


4.857 


4. 


3.4 


3 




45.89 


45.89 


44.92 


30.22 


21.22 


16 


^5 


656.3 


654.6 


624.4 


342.2 


196.4 


125 


-^6 


12,653 


12,568 


11,666 


5200 


2429 


1296 



Table 2: Moments of S field 





Exact 


PPZA 


PZA 


ZA 


LP 


FF 




-3.714 


-3.714 


-3.714 


-2 


-0.8 


0. 




27.41 


27.41 


24.49 


8. 


1.46 


0. 




-308.4 


-301.5 


-240.8 


-48.9 


-4.19 


0. 




4694 


4450 


3180 


404.4 


16.35 


0. 



Table 3: Moments of 6 field 



18 



References 

Bagla, J.S., & Padmanabhan, T., 1994. MNRAS, 266, 227. 



[1] 
[2] 
[3] 
[4] 
[5 



Bernardeau, F., 1992. ApJ, 392, 1. 

Bernardeau, F., 1993. CITA preprint 93/14; ApJ, in press. 



Bernardeau, F., & Schaeffer, R., 1992. ALA, 255, 1. 



Bernardeau, F., Singh, T.P., Banerjee B.,& Chitre S.M., 1993. Preprint TIFR, \astro- 
ph/9311053i ; MNRAS, submitted. 



[6] Bouchet, F.R, Juszkiewicz, R., Colombi, S. & Pellat, R.,1992. ApJL, 394, L5. 
[7] Brainerd, T.G., Scherer, R.J., & Villumsen, J.V., 1993. ApJ, 418, 570. 

Buchert, T., 1992. MNRAS, 254, 729. 
[9] Coles, P., & Jones, B., 1991. MNRAS, 248, 1. 
[10] Fry, J.N., 1984. ApJ, 279, 499. 
[11] Gramman, M., 1993. ApJL, 405, L47. 
[12] Grinstein, B., & Wise, M.B., 1987. ApJ, 320, 448. 
[13] Hamilton, A., 1998. ApJL, 331, L59. 

[14] Juszkiewicz, R., Weiberg, D.H., Amsterdamsky, P., Chodorovsky, M., & Bouchet, F.R., 
1993. IAS preprint (lASSNS-AST 93/50). 

[15] Kofman, L., 1991. In: Primordial Nucleosynthesis and Evolution of Early Universe, eds. 
K.Sato & J. Audouze Dordrecht: Kluwer, p. 495. 

19 



[16] Kofman, L., Bertschinger, E., Gelb, J.M., Nusser, A., & Dekel, A., 1994. ApJ, 420, 44. 
[17] Lachieze-Rey, M., 1993. ApJ 408, 403. 

[18] Matarrese, S., Lucchin, F., Moscardini, L., & Saez, D., 1992. MNRAS, 259, 437. 



[19] Moutarde, F., Alimi, J.M., Bouchet, F.R., Pellat, R. & Ramani, A., 1991. ApJ, 382, 
377. 



[20] Munshi, D., & Starobinsky, A.A., 1993. Preprint lUCAA, \astro-ph/ 93110501 ; ApJ, in 
press. 

[21] Padmanabhan, T. & Subramanian, K., 1993. Ap.J, 410, 482. 

[22] Peebles, P.J.E., 1980. The Large-Scale Structure of the Universe, Princeton, Princeton 
University Press. 

[23] Sathyaprakash, B.S., Munshi, D., Sahni, V., Pogosyan, D. & Melott, A.L., 1994, in 
preparation. 

[24] Shandarin, S.F., & Zeldovich, Ya. B., 1989. Rev. Mod. Phys., 61, 185. 

[25] White, S.D.M., 1979. MNRAS, 186, 145. 

[26] Zeldovich, Ya.B., 1970. Astron.Astroph., 5, 84. 



[27] Zeldovich, Ya.B. and Novikov, I.D., 1983, The Structure and Evolution of the Universe 
(University of Chicago, Chicago/London). 



20 



Figure Captions 

Fig. 1 : The density in the top-hat spherical collapse model as estimated in the different 
approximations {y{r) = 1 + 5) is plotted against the exact solution. The exact solution is 
labelled by 1; PPZA by 2; PZA by 3; ZA by 4; LP by 5; FF by 6; and the linear solution by 
7. 
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This figure "figl-l.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/9402065vl 



